home *** CD-ROM | disk | FTP | other *** search
/ ADA Programming Guide / ADA Programming Guide.iso / ada_gnu / adainc / a-ngcoty.adb < prev    next >
Text File  |  1996-01-30  |  14KB  |  535 lines

  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT RUNTIME COMPONENTS                          --
  4. --                                                                          --
  5. --   A D A . N U M E R I C S . G E N E R I C _ C O M P L E X _ T Y P E S    --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.4 $                              --
  10. --                                                                          --
  11. --           Copyright (c) 1992,1993,1994 NYU, All Rights Reserved          --
  12. --                                                                          --
  13. -- The GNAT library is free software; you can redistribute it and/or modify --
  14. -- it under terms of the GNU Library General Public License as published by --
  15. -- the Free Software  Foundation; either version 2, or (at your option) any --
  16. -- later version.  The GNAT library is distributed in the hope that it will --
  17. -- be useful, but WITHOUT ANY WARRANTY;  without even  the implied warranty --
  18. -- of MERCHANTABILITY  or  FITNESS FOR  A PARTICULAR PURPOSE.  See the  GNU --
  19. -- Library  General  Public  License for  more  details.  You  should  have --
  20. -- received  a copy of the GNU  Library  General Public License  along with --
  21. -- the GNAT library;  see the file  COPYING.LIB.  If not, write to the Free --
  22. -- Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.        --
  23. --                                                                          --
  24. ------------------------------------------------------------------------------
  25.  
  26. with Ada.Numerics.Aux; use Ada.Numerics.Aux;
  27. package body Ada.Numerics.Generic_Complex_Types is
  28.  
  29.    subtype R is Real'Base;
  30.  
  31.    ---------
  32.    -- "+" --
  33.    ---------
  34.  
  35.    function "+" (Right : Complex) return Complex is
  36.    begin
  37.       return Right;
  38.    end "+";
  39.  
  40.    function "+" (Left, Right : Complex) return Complex is
  41.    begin
  42.       return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
  43.    end "+";
  44.  
  45.    function "+" (Right : Imaginary) return Imaginary is
  46.    begin
  47.       return Right;
  48.    end "+";
  49.  
  50.    function "+" (Left, Right : Imaginary) return Imaginary is
  51.    begin
  52.       return Imaginary (R (Left) + R (Right));
  53.    end "+";
  54.  
  55.    function "+" (Left : Complex; Right : Real'Base) return Complex is
  56.    begin
  57.       return Complex'(Left.Re + Right, Left.Im);
  58.    end "+";
  59.  
  60.    function "+" (Left : Real'Base; Right : Complex) return Complex is
  61.    begin
  62.       return Complex'(Left + Right.Re, Right.Im);
  63.    end "+";
  64.  
  65.    function "+" (Left : Complex; Right : Imaginary) return Complex is
  66.    begin
  67.       return Complex'(Left.Re, Left.Im + R (Right));
  68.    end "+";
  69.  
  70.    function "+" (Left : Imaginary; Right : Complex) return Complex is
  71.    begin
  72.       return Complex'(R (Left) + Right.Re, Right.Im);
  73.    end "+";
  74.  
  75.    function "+" (Left : Imaginary; Right : Real'Base) return Complex is
  76.    begin
  77.       return Complex'(Right, R (Left));
  78.    end "+";
  79.  
  80.    function "+" (Left : Real'Base; Right : Imaginary) return Complex is
  81.    begin
  82.       return Complex'(Left, R (Right));
  83.    end "+";
  84.  
  85.    ---------
  86.    -- "-" --
  87.    ---------
  88.  
  89.    function "-" (Right : Complex) return Complex is
  90.    begin
  91.       return (-Right.Re, -Right.Im);
  92.    end "-";
  93.  
  94.    function "-" (Left, Right : Complex) return Complex is
  95.    begin
  96.       return (Left.Re - Right.Re, Left.Im - Right.Im);
  97.    end "-";
  98.  
  99.    function "-" (Right : Imaginary) return Imaginary is
  100.    begin
  101.       return Imaginary (-R (Right));
  102.    end "-";
  103.  
  104.    function "-" (Left, Right : Imaginary) return Imaginary is
  105.    begin
  106.       return Imaginary (R (Left) - R (Right));
  107.    end "-";
  108.  
  109.    function "-" (Left : Complex; Right : Real'Base) return Complex is
  110.    begin
  111.       return Complex'(Left.Re - Right, Left.Im);
  112.    end "-";
  113.  
  114.    function "-" (Left : Real'Base; Right : Complex) return Complex is
  115.    begin
  116.       return Complex'(Left - Right.Re, -Right.Im);
  117.    end "-";
  118.  
  119.    function "-" (Left : Complex; Right : Imaginary) return Complex is
  120.    begin
  121.       return Complex'(Left.Re, Left.Im - R (Right));
  122.    end "-";
  123.  
  124.    function "-" (Left : Imaginary; Right : Complex) return Complex is
  125.    begin
  126.       return Complex'(R (Left) - Right.Re, -Right.Im);
  127.    end "-";
  128.  
  129.    function "-" (Left : Imaginary; Right : Real'Base) return Complex is
  130.    begin
  131.       return Complex'(-Right, R (Left));
  132.    end "-";
  133.  
  134.    function "-" (Left : Real'Base; Right : Imaginary) return Complex is
  135.    begin
  136.       return Complex'(Left, -R (Right));
  137.    end "-";
  138.  
  139.    ---------
  140.    -- "*" --
  141.    ---------
  142.  
  143.    function "*" (Left, Right : Complex) return Complex is
  144.    begin
  145.       return  (Re => Left.Re * Right.Re - Left.Im * Right.Im,
  146.                Im => Left.Re * Right.Im + Left.Im * Right.Re);
  147.    end "*";
  148.  
  149.    function "*" (Left, Right : Imaginary) return Real'Base is
  150.    begin
  151.       return -R (Left) * R (Right);
  152.    end "*";
  153.  
  154.    function "*" (Left : Complex; Right : Real'Base) return Complex is
  155.    begin
  156.       return Complex'(Left.Re * Right, Left.Im * Right);
  157.    end "*";
  158.  
  159.    function "*" (Left : Real'Base; Right : Complex) return Complex is
  160.    begin
  161.       return (Left * Right.Re, Left * Right.Im);
  162.    end "*";
  163.  
  164.    function "*" (Left : Complex; Right : Imaginary) return Complex is
  165.    begin
  166.       return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
  167.    end "*";
  168.  
  169.    function "*" (Left : Imaginary; Right : Complex) return Complex is
  170.    begin
  171.       return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
  172.    end "*";
  173.  
  174.    function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
  175.    begin
  176.       return Left * Imaginary (Right);
  177.    end "*";
  178.  
  179.    function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
  180.    begin
  181.       return Imaginary (Left * R (Right));
  182.    end "*";
  183.  
  184.    ---------
  185.    -- "/" --
  186.    ---------
  187.  
  188.    function "/" (Left, Right : Complex) return Complex is
  189.       a : constant R := Left.Re;
  190.       b : constant R := Left.Im;
  191.       c : constant R := Right.Re;
  192.       d : constant R := Right.Im;
  193.  
  194.    begin
  195.       return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
  196.                       Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
  197.    end "/";
  198.  
  199.    function "/" (Left, Right : Imaginary) return Real'Base is
  200.    begin
  201.       return R (Left) / R (Right);
  202.    end "/";
  203.  
  204.    function "/" (Left : Complex; Right : Real'Base) return Complex is
  205.    begin
  206.       return Complex'(Left.Re / Right, Left.Im / Right);
  207.    end "/";
  208.  
  209.    function "/" (Left : Real'Base; Right : Complex) return Complex is
  210.       a : constant R := Left;
  211.       c : constant R := Right.Re;
  212.       d : constant R := Right.Im;
  213.    begin
  214.       return Complex'(Re =>  (a * c) / (c ** 2 + d ** 2),
  215.                       Im => -(a * d) / (c ** 2 + d ** 2));
  216.    end "/";
  217.  
  218.    function "/" (Left : Complex; Right : Imaginary) return Complex is
  219.       a : constant R := Left.Re;
  220.       b : constant R := Left.Im;
  221.       d : constant R := R (Right);
  222.  
  223.    begin
  224.       return (b / d,  -a / d);
  225.    end "/";
  226.  
  227.    function "/" (Left : Imaginary; Right : Complex) return Complex is
  228.       b : constant R := R (Left);
  229.       c : constant R := Right.Re;
  230.       d : constant R := Right.Im;
  231.  
  232.    begin
  233.       return (Re => -b * d / (c ** 2 + d ** 2),
  234.               Im => b * c / (c ** 2 + d ** 2));
  235.    end "/";
  236.  
  237.    function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
  238.    begin
  239.       return Imaginary (R (Left) / Right);
  240.    end "/";
  241.  
  242.    function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
  243.    begin
  244.       return Imaginary (-Left / R (Right));
  245.    end "/";
  246.  
  247.    ----------
  248.    -- "**" --
  249.    ----------
  250.  
  251.    function "**" (Left : Complex; Right : Integer) return Complex is
  252.       Result : Complex := (1.0, 0.0);
  253.       Factor : Complex := Left;
  254.       Exp    : Natural := Right;
  255.  
  256.    begin
  257.       --  We use the standard logarithmic approach, Exp gets shifted right
  258.       --  testing successive low order bits and Factor is the value of the
  259.       --  base raised to the next power of 2. For positive exponents we
  260.       --  multiply the result by this factor, for negative exponents, we
  261.       --  divide by this factor.
  262.  
  263.       if Exp >= 0 then
  264.  
  265.          --  For a positive exponent, if we get a constraint error during
  266.          --  this loop, it is an overflow, and the constraint error will
  267.          --  simply be passed on to the caller.
  268.  
  269.          while Exp /= 0 loop
  270.             if Exp rem 2 /= 0 then
  271.                Result := Result * Factor;
  272.             end if;
  273.  
  274.             Factor := Factor * Factor;
  275.             Exp := Exp / 2;
  276.          end loop;
  277.  
  278.          return Result;
  279.  
  280.       else -- Exp < 0 then
  281.  
  282.          --  For the negative exponent case, a constraint error during this
  283.          --  calculation happens if Factor gets too large, and the proper
  284.          --  response is to return 0.0, since what we essentially have is
  285.          --  1.0 / infinity, and the closest model number will be zero.
  286.  
  287.          begin
  288.  
  289.             while Exp /= 0 loop
  290.                if Exp rem 2 /= 0 then
  291.                   Result := Result * Factor;
  292.                end if;
  293.  
  294.                Factor := Factor * Factor;
  295.                Exp := Exp / 2;
  296.             end loop;
  297.  
  298.             return R ' (1.0) / Result;
  299.  
  300.          exception
  301.  
  302.             when Constraint_Error =>
  303.                return (0.0, 0.0);
  304.          end;
  305.       end if;
  306.    end "**";
  307.  
  308.    function "**" (Left : Imaginary; Right : Integer) return Complex is
  309.       M : R := R (Left) ** Right;
  310.    begin
  311.       case Right mod 4 is
  312.          when 0 => return (M,   0.0);
  313.          when 1 => return (0.0, M);
  314.          when 2 => return (-M,  0.0);
  315.          when 3 => return (0.0, -M);
  316.          when others => raise Program_Error;
  317.       end case;
  318.    end "**";
  319.  
  320.    ---------
  321.    -- "<" --
  322.    ---------
  323.  
  324.    function "<" (Left, Right : Imaginary) return Boolean is
  325.    begin
  326.       return R (Left) < R (Right);
  327.    end "<";
  328.  
  329.    ----------
  330.    -- "<=" --
  331.    ----------
  332.  
  333.    function "<=" (Left, Right : Imaginary) return Boolean is
  334.    begin
  335.       return R (Left) <= R (Right);
  336.    end "<=";
  337.  
  338.    ---------
  339.    -- ">" --
  340.    ---------
  341.  
  342.    function ">" (Left, Right : Imaginary) return Boolean is
  343.    begin
  344.       return R (Left) > R (Right);
  345.    end ">";
  346.  
  347.    ----------
  348.    -- ">=" --
  349.    ----------
  350.  
  351.    function ">=" (Left, Right : Imaginary) return Boolean is
  352.    begin
  353.       return R (Left) >= R (Right);
  354.    end ">=";
  355.  
  356.    -----------
  357.    -- "abs" --
  358.    -----------
  359.  
  360.    function "abs" (Right : Imaginary) return Real'Base is
  361.    begin
  362.       return R (Right);
  363.    end "abs";
  364.  
  365.    --------------
  366.    -- Argument --
  367.    --------------
  368.  
  369.    function Argument (X : Complex) return Real'Base is
  370.       a : constant R := X.Re;
  371.       b : constant R := X.Im;
  372.  
  373.    begin
  374.       if b = 0.0 then
  375.          if a >= 0.0 then
  376.             return 0.0;
  377.          else
  378.             return Pi;
  379.          end if;
  380.       else
  381.          return R (Atan (Double (a / b)));
  382.       end if;
  383.  
  384.    exception
  385.       when Constraint_Error =>
  386.          if a > 0.0 then
  387.             return 0.0;
  388.          else
  389.             return Pi;
  390.          end if;
  391.    end Argument;
  392.  
  393.    function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
  394.    begin
  395.       if Cycle > 0.0 then
  396.          return Argument (X) * Cycle / (2.0 * Pi);
  397.       else
  398.          raise Constraint_Error;
  399.       end if;
  400.    end Argument;
  401.  
  402.    ----------------------------
  403.    -- Compose_From_Cartesian --
  404.    ----------------------------
  405.  
  406.    function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
  407.    begin
  408.       return (Re, Im);
  409.    end Compose_From_Cartesian;
  410.  
  411.    function Compose_From_Cartesian (Re : Real'Base) return Complex is
  412.    begin
  413.       return (Re, 0.0);
  414.    end Compose_From_Cartesian;
  415.  
  416.    function Compose_From_Cartesian (Im : Imaginary) return Complex is
  417.    begin
  418.       return (0.0, R (Im));
  419.    end Compose_From_Cartesian;
  420.  
  421.    ------------------------
  422.    -- Compose_From_Polar --
  423.    ------------------------
  424.  
  425.    function Compose_From_Polar (
  426.      Modulus, Argument : Real'Base)
  427.      return Complex
  428.    is
  429.    begin
  430.       if Modulus = 0.0 then
  431.          return (0.0, 0.0);
  432.       else
  433.          return (Modulus * R (Cos (Double (Argument))),
  434.                  Modulus * R (Sin (Double (Argument))));
  435.       end if;
  436.    end Compose_From_Polar;
  437.  
  438.    function Compose_From_Polar (
  439.      Modulus, Argument, Cycle : Real'Base)
  440.      return Complex
  441.    is
  442.       Arg : Real'Base;
  443.  
  444.    begin
  445.       if Modulus = 0.0 then
  446.          return (0.0, 0.0);
  447.  
  448.       elsif Cycle > 0.0 then
  449.          if Argument = 0.0 then
  450.             return (Modulus, 0.0);
  451.  
  452.          elsif Argument = Cycle / 4.0 then
  453.             return (0.0, Modulus);
  454.  
  455.          elsif Argument = Cycle / 2.0 then
  456.             return (-Modulus, 0.0);
  457.  
  458.          elsif Argument = 3.0 * Cycle / 4.0 then
  459.             return (0.0, -Modulus);
  460.          else
  461.             Arg := 2.0 * Pi * Argument / Cycle;
  462.             return (Modulus * R (Cos (Double (Arg))),
  463.                     Modulus * R (Sin (Double (Arg))));
  464.          end if;
  465.       else
  466.          raise Constraint_Error;
  467.       end if;
  468.    end Compose_From_Polar;
  469.  
  470.    ---------------
  471.    -- Conjugate --
  472.    ---------------
  473.  
  474.    function Conjugate (X : Complex) return Complex is
  475.    begin
  476.       return Complex'(X.Re, -X.Im);
  477.    end Conjugate;
  478.  
  479.    --------
  480.    -- Im --
  481.    --------
  482.  
  483.    function Im (X : Complex) return Real'Base is
  484.    begin
  485.       return X.Im;
  486.    end Im;
  487.  
  488.    function Im (X : Imaginary) return Real'Base is
  489.    begin
  490.       return R (X);
  491.    end Im;
  492.  
  493.    -------------
  494.    -- Modulus --
  495.    -------------
  496.  
  497.    function Modulus (X : Complex) return Real'Base is
  498.    begin
  499.       return R (Sqrt (Double (X.Re ** 2 + X.Im ** 2)));
  500.    end Modulus;
  501.  
  502.    --------
  503.    -- Re --
  504.    --------
  505.  
  506.    function Re (X : Complex) return Real'Base is
  507.    begin
  508.       return X.Re;
  509.    end Re;
  510.  
  511.    ------------
  512.    -- Set_Im --
  513.    ------------
  514.  
  515.    procedure Set_Im (X : in out Complex; Im : in Real'Base) is
  516.    begin
  517.       X.Im := Im;
  518.    end Set_Im;
  519.  
  520.    procedure Set_Im (X : out Imaginary; Im : in Real'Base) is
  521.    begin
  522.       X := Imaginary (Im);
  523.    end Set_Im;
  524.  
  525.    ------------
  526.    -- Set_Re --
  527.    ------------
  528.  
  529.    procedure Set_Re (X : in out Complex; Re : in Real'Base) is
  530.    begin
  531.       X.Re := Re;
  532.    end Set_Re;
  533.  
  534. end Ada.Numerics.Generic_Complex_Types;
  535.